Aplicaciones de Agro Data Science

Análisis de mortalidad de insectos y predicción de rendimiento

Christian R.A. Vásquez Velasco

InkaStats Academy

2024-09-15

Agro Data Science

La agricultura moderna

Agro Data Science

ADS combina el conocimiento agronómico con las tecnologías de la información para resolver problemas en la agricultura.
Abarca temas desde la recopilación y análisis de datos hasta la visualización y comunicación de resultados.
Los científicos de datos agrícolas trabajan con agricultores, empresas agroindustriales, universidades y organizaciones gubernamentales para mejorar la eficiencia y la sostenibilidad de la agricultura.

Problemas que resuelven el ADS

Gestión del riesgo climático
Determinación de dosis y momento de fertilización
Análisis y predicción de rendimiento
Estrés hídrico
Monitoreo de campos
Diagnóstico de plagas y enfermedades

Importancia del acceso a datos precisos y confiables: Estaciones meteorológicas, sensores de campo, drones y satélites.
Herramientas y técnicas de análisis de datos: Visualización de datos, Clasificación y Regresión, Series temporales, Aprendizaje automático, Inteligencia artificial.

Desafios

Falta de acceso a datos de calidad.
Cuestiones éticas y de privacidad.

Oportunidades

Mejorar la eficiencia y la sostenibilidad de la agricultura.
Tomar decisiones informadas y basadas en datos.

R

R es un lenguaje de programación utilizado en el análisis estadístico y la visualización de datos.
Fue creado en 1993 por Ross Ihaka y Robert Gentleman en la Universidad de Auckland, Nueva Zelanda.
R es un software libre y de código abierto.
R cuenta con una gran cantidad de paquetes que se utilizan para realizar tareas específicas en el análisis de datos.
R tiene una sintaxis fácil de aprender y es utilizado tanto por científicos de datos como por estadísticos, ingenieros y profesionales de negocios.

Tidyverse

Tidyverse es un conjunto de paquetes de software de análisis de datos en el lenguaje de programación R, que se utilizan para manipular, visualizar y modelar datos.
Estos paquetes comparten una filosofía común y una sintaxis consistente, lo que hace que sea fácil trabajar con ellos en conjunto.
Los paquetes que forman parte de Tidyverse incluyen dplyr, ggplot2, tidyr, readr, purrr, tibble, stringr, forcats, entre otros.

Tidymodels

Tidymodels es un conjunto de paquetes de software de modelado estadístico en el lenguaje de programación R que se utilizan para construir, ajustar, validar y analizar modelos predictivos.
El objetivo principal de Tidymodels es proporcionar herramientas para ayudar a los usuarios a realizar modelado estadístico de manera más eficiente y efectiva, utilizando técnicas modernas de modelado y validación..
Los paquetes que forman parte de Tidymodels incluyen recipes, parsnip, yardstick, dials, infer, entre otros.

Análisis de mortalidad de insectos

Entendiendo el problema

Las áreas de sanidad cuentan con data histórica.

Deben transformarla en información útil.

¿?

Predecir la mortalidad que se obtendrá en las aplicaciones futuras.

Para ilustración se utilizará el inventario de Eficacia de plaguicidas.

Data summary
Name	data
Number of rows	929
Number of columns	6
_______________________
Column type frequency:
factor	3
numeric	3
________________________
Group variables	None

Variable type: factor

skim_variable	complete_rate	ordered	n_unique	top_counts
InsectStage	1	FALSE	134	133: 41, 134: 41, 62: 39, 61: 32
Method	1	FALSE	7	5: 707, 1: 133, 4: 75, 6: 6
Product	1	FALSE	101	1: 145, 94: 38, 64: 26, 73: 26

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p25	p50	p75	p100	hist
Pre_Mean_Count	0	1.00	83.41	177.35	6.41	24.23	85.67	2416	▇▁▁▁▁
Post_Mean_Count	5	0.99	32.05	123.41	0.22	2.33	18.19	1929	▇▁▁▁▁
Mortality_Mean_Cor	0	1.00	0.63	0.41	0.15	0.84	0.98	1	▃▁▁▂▇

Método: método de aplicación
- 1 Control - 2 Insecticidas en gránulos - 3 Tratamiento de semillas - 4 Aplicación al suelo - 5 Pulverización - 6 Niebla de vapor - 7 Aplicación superficial con agua

# A tibble: 10 × 6
   InsectStage Method Product Pre_Mean_Count Post_Mean_Count Mortality_Mean_Cor
   <fct>       <fct>  <fct>            <dbl>           <dbl>              <dbl>
 1 1           1      1                87.8           150.                0    
 2 1           3      57                1.67            3                 0    
 3 1           3      71                2.67            2.67              0    
 4 1           4      31               45               6.33              0.859
 5 1           4      49               45               5.67              0.874
 6 1           4      57               90.6            48.6               0.464
 7 1           4      78               39.3            16.3               0.585
 8 1           5      58               74.3            46.7               0.372
 9 2           1      1                13.3            53.7               0    
10 2           4      57                8              11                 0

Partición inicial

split <- initial_split(data, prop = 0.70, strata = Mortality_Mean_Cor)

Partición inicial

split <- initial_split(data, prop = 0.70, strata = Mortality_Mean_Cor)
data_train <- split %>% training()

tibble [649 × 6] (S3: tbl_df/tbl/data.frame)
 $ InsectStage       : Factor w/ 140 levels "1","2","3","4",..: 1 2 2 2 3 3 4 4 5 7 ...
 $ Method            : Factor w/ 7 levels "1","2","3","4",..: 1 1 4 5 1 3 1 4 1 1 ...
 $ Product           : Factor w/ 101 levels "1","2","3","4",..: 1 1 57 58 1 57 1 58 1 1 ...
 $ Pre_Mean_Count    : num [1:649] 87.8 13.3 8 14.7 117.3 ...
 $ Post_Mean_Count   : num [1:649] 150.5 53.7 11 18.3 167.9 ...
 $ Mortality_Mean_Cor: num [1:649] 0 0 0 0 0 ...

Partición inicial

split <- initial_split(data, prop = 0.70, strata = Mortality_Mean_Cor)
data_train <- split %>% training()
data_test <- split %>% testing()

tibble [280 × 6] (S3: tbl_df/tbl/data.frame)
 $ InsectStage       : Factor w/ 140 levels "1","2","3","4",..: 1 1 1 1 1 3 4 6 6 6 ...
 $ Method            : Factor w/ 7 levels "1","2","3","4",..: 3 3 4 4 4 3 4 1 5 5 ...
 $ Product           : Factor w/ 101 levels "1","2","3","4",..: 57 71 49 57 78 71 57 1 14 69 ...
 $ Pre_Mean_Count    : num [1:280] 1.67 2.67 45 90.6 39.33 ...
 $ Post_Mean_Count   : num [1:280] 3 2.67 5.67 48.6 16.33 ...
 $ Mortality_Mean_Cor: num [1:280] 0 0 0.874 0.464 0.585 ...

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train)

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`))

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("Mortality_Mean_Cor")),
                 fn = ~log(.))

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("Mortality_Mean_Cor")),
                 fn = ~log(.))

Receta de preprocesamiento

data_rec <- recipe(formula = `Mortality_Mean_Cor` ~ .,
                   x = data_train) %>%
  step_select(-c(`Post_Mean_Count`)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("Mortality_Mean_Cor")),
                 fn = ~log(.)) %>%
  step_dummy(Method, one_hot = TRUE) %>%
  step_dummy(Product, one_hot = TRUE) %>%
  step_dummy(InsectStage, one_hot = TRUE)

Receta de preprocesamiento

data_prep <- prep(data_rec)

Receta de preprocesamiento

data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)

Receta de preprocesamiento

data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)
data_test_prep  <- bake(data_prep, new_data = data_test)

Receta de preprocesamiento

data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)
data_test_prep  <- bake(data_prep, new_data = data_test)
head(data_train_prep)

# A tibble: 6 × 250
  Pre_Mean_Count Mortality_Mean_Cor Method_X1 Method_X2 Method_X3 Method_X4
           <dbl>              <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
1           4.48                  0         1         0         0         0
2           2.59                  0         1         0         0         0
3           2.08                  0         0         0         0         1
4           2.69                  0         0         0         0         0
5           4.76                  0         1         0         0         0
6           2.46                  0         0         0         1         0
# ℹ 244 more variables: Method_X5 <dbl>, Method_X6 <dbl>, Method_X7 <dbl>,
#   Product_X1 <dbl>, Product_X2 <dbl>, Product_X3 <dbl>, Product_X4 <dbl>,
#   Product_X5 <dbl>, Product_X6 <dbl>, Product_X7 <dbl>, Product_X8 <dbl>,
#   Product_X9 <dbl>, Product_X10 <dbl>, Product_X11 <dbl>, Product_X12 <dbl>,
#   Product_X13 <dbl>, Product_X14 <dbl>, Product_X15 <dbl>, Product_X16 <dbl>,
#   Product_X17 <dbl>, Product_X18 <dbl>, Product_X19 <dbl>, Product_X20 <dbl>,
#   Product_X21 <dbl>, Product_X22 <dbl>, Product_X23 <dbl>, …

Random Forest

Definir especificaciones del modelo

# Inicializar un objeto de random forest
rf_model_spec <- rand_forest(
  trees = tune(), 
  min_n = tune(), 
  mtry = tune()
) %>% 
  # Establecer el modelo de motor
  set_engine('ranger',                         
             importance = "permutation") %>% 
  # Establecer el modo de modelo
  set_mode('regression')

Definir fórmula del modelo

formula_rf = Mortality_Mean_Cor ~ .

Definir flujo de trabajo

rf_wf <- workflow() %>%
  add_recipe(data_rec) %>%
  add_model(rf_model_spec,
            formula = formula_rf) 
#obtener todos los argumentos ajustables posibles en el flujo de trabajo
tune_args(rf_wf)

# A tibble: 3 × 6
  name  tunable id    source     component   component_id
  <chr> <lgl>   <chr> <chr>      <chr>       <chr>       
1 mtry  TRUE    mtry  model_spec rand_forest <NA>        
2 trees TRUE    trees model_spec rand_forest <NA>        
3 min_n TRUE    min_n model_spec rand_forest <NA>

Extraer el conjunto de hiperparámetros del modelo

rf_params <- hardhat::extract_parameter_set_dials(rf_wf) %>% 
  update(trees = trees(range = c(50L, 150L)))
rf_params <- rf_params %>% 
  dials::finalize(rf_params)
rf_params

Collection of 3 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[+]
      trees trees nparam[+]
      min_n min_n nparam[+]

Definir métricas del performance predictivo

multi_met <- yardstick::metric_set(
yardstick::rmse,
yardstick::rsq,
yardstick::mape,
yardstick::mae)

Definir esquema de validación cruzada

set.seed(234)
data_folds <- rsample::vfold_cv(data_train,
                                 v = 20)

Búsqueda del mejor conjunto de Hiperparámtros mediante validación cruzada

set.seed(2020)
tictoc::tic()
rf_tune <-
  rf_wf %>%
  tune_bayes(
    resamples = data_folds,
    param_info = rf_params,
    # por defecto es 5 (número aleatorio de combinaciones (puntos de grilla) de hiperparámetros)
    initial = 5, 
    iter = 50,
    metrics = multi_met,
    control = control_bayes(
      # El corte entero para el número de iteraciones sin mejores resultados.
      no_improve = 10,
      extract = identity,
      save_pred = TRUE,
      verbose = FALSE
      )
  )
tictoc::toc()

320 sec elapsed

Métricas del rendimiento predictivo en fase de validación cruzada

Seleccionar el mejor conjunto de hiperparámetros

best_rf_model <- select_best(rf_tune, "rmse")
final_rf <- finalize_workflow(rf_wf,
                              best_rf_model)
final_rf

══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_select()
• step_mutate_at()
• step_dummy()
• step_dummy()
• step_dummy()

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (regression)

Main Arguments:
  mtry = 6
  trees = 76
  min_n = 39

Engine-Specific Arguments:
  importance = permutation

Computational engine: ranger

Métricas del rendimiento predictivo para el mejor conjunto de hiperparámetros

# Métricas promedio de todas las particiones
metrics_rf <- 
  rf_tune %>% 
  collect_metrics(summarize = TRUE) %>%
  dplyr::mutate(modelo = "rf") %>%
  dplyr::select(-c(.estimator,.config)) %>%
  dplyr::filter(mtry %in% best_rf_model$mtry,
                min_n %in% best_rf_model$min_n,
                trees %in% best_rf_model$trees)
metrics_rf

mtry	trees	min_n	.metric	mean	n	std_err	.iter	modelo
6	76	39	mae	0.2134	20	0.004847	11	rf
6	76	39	mape	Inf	20	NaN	11	rf
6	76	39	rmse	0.2688	20	0.006324	11	rf
6	76	39	rsq	0.5735	20	0.024704	11	rf

Entrenamiento del modelo final

rf_model <- final_rf %>%
  fit(data_train) %>%
  extract_fit_parsnip()
rf_model

parsnip model object

Ranger result

Call:
 ranger::ranger(x = maybe_data_frame(x), y = y, mtry = min_cols(~6L,      x), num.trees = ~76L, min.node.size = min_rows(~39L, x),      importance = ~"permutation", num.threads = 1, verbose = FALSE,      seed = sample.int(10^5, 1)) 

Type:                             Regression 
Number of trees:                  76 
Sample size:                      649 
Number of independent variables:  249 
Mtry:                             6 
Target node size:                 39 
Variable importance mode:         permutation 
Splitrule:                        variance 
OOB prediction error (MSE):       0.07622 
R squared (OOB):                  0.5352

Importancia de variables

Testeo interno

Resumen de la predicción

predicciones <- rf_model %>%
  predict(
    new_data = data_test_prep,
    type = "numeric"
  )
  
predicciones_rf <- predicciones %>% 
  bind_cols(data_test_prep %>% dplyr::select(Mortality_Mean_Cor)) %>%
  dplyr::mutate(modelo = "rf")

summary(predicciones_rf)

     .pred       Mortality_Mean_Cor    modelo         
 Min.   :0.127   Min.   :0.000      Length:280        
 1st Qu.:0.530   1st Qu.:0.135      Class :character  
 Median :0.719   Median :0.847      Mode  :character  
 Mean   :0.635   Mean   :0.623                        
 3rd Qu.:0.809   3rd Qu.:0.983                        
 Max.   :0.919   Max.   :1.000

Testeo interno

Resumen de las métricas de predicción

error_test_rf  <- multi_met(
  data     = predicciones_rf,
  truth    = Mortality_Mean_Cor,
  estimate = .pred,
  na_rm    = TRUE
) %>%
  dplyr::mutate(
    modelo = "rf"
  )
error_test_rf

# A tibble: 4 × 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard       0.275 rf    
2 rsq     standard       0.591 rf    
3 mape    standard     Inf     rf    
4 mae     standard       0.225 rf

Testeo interno

Predicción de rendimiento

Fórmulas para la proyección del rendimiento

Debemos entender que de forma empírica, las estimaciones de cosecha se realizan en muchos cultivos teniendo en cuenta una fórmula matemática como la siguiente:

\[Producción = A * DS * NOP * PO\]

Dónde:

A = Área

DS = Densidad de siembra

NOP = Número de órganos por planta

PO = Peso de órganos

Entendiéndose con esta fórmula que para estimar la producción por campaña de cualquier cultivo, se necesita información de la extensión de siembra, número de plantas por hectárea, número de órganos a cosechar y el peso promedio de estos órganos.

Para fines logísticos, hacer una estimación de producción a nivel de campaña no ayuda mucho, ya que se necesita saber la producción en distintos puntos de corte o periodos de tiempo.

Para esto último fin, podemos hacer una modificación a la ecuación, obteniendo lo siguiente:

\[Producción_{(t)} = A * DS * NOP_{(t)} * PO_{(t)} * TM_{(t)}\]

Dónde:

A = Área

DS = Densidad de siembra

NOP = Número de órganos por planta

PO = Peso de órganos

TM = Tasa de maduración

Esta última ecuación añade un modelamiento temporal, que está en función de la tasa de maduración y crecimiento de los órganos a cosechar.

Entonces, para realizar Crop Yield Forecast, debemos entender que existe en el tiempo uno o más factores que condicionarán la producción en distintos puntos transversales y por ello debemos tomar datos horizontales o de series de tiempo en el cultivo, pues serán fundamentales para realizar las proyecciones. Es así, que debemos imaginar cualquier algoritmo como una función similar a esta última ecuación matemática y debemos lograr que el modelo por crear, a través del aprendizaje estadístico, cumpla con la lógica necesaria para simular este comportamiento matemático (que fue pensado para simular el comportamiento real del cultivo).

Wild Blueberry Yield Prediction

El artículo “Wild blueberry yield prediction using a combination of computer simulation and machine learning algorithms”, usa datos recolectados por 30 años en Maine (EEUU) y publicados por Wild Blueberry Pollination Model, para un modelo de simulación espacial. El modelamiento predictivo en arándanos requiere de datos con fuerte influyencia de factores espaciales (geográficos), y adicionalmente en este caso, de las plantas, especies de abejas y meteorología. Estos datos simulados son usados por investigadores para probar la calidad predictiva de sus algoritmos de aprendizaje estadístico a datos reales y datos generados por simulación. Emplearemos este conjunto de datos para entrenar modelos en el enfoque tidy.

Descriptive Statistics  

                    andrena   AverageOfLowerTRange   AverageOfUpperTRange   AverageRainingDays
----------------- --------- ---------------------- ---------------------- --------------------
             Mean      0.47                  48.61                  68.72                 0.32
          Std.Dev      0.16                   5.42                   7.68                 0.17
              Min      0.00                  41.20                  58.20                 0.06
               Q1      0.38                  45.80                  64.70                 0.10
           Median      0.50                  50.80                  71.90                 0.26
               Q3      0.63                  50.80                  71.90                 0.39
              Max      0.75                  55.90                  79.00                 0.56
              MAD      0.18                   7.41                  10.53                 0.24
              IQR      0.25                   5.00                   7.20                 0.29
               CV      0.34                   0.11                   0.11                 0.54
         Skewness      0.19                  -0.01                  -0.02                 0.07
      SE.Skewness      0.09                   0.09                   0.09                 0.09
         Kurtosis     -0.67                  -1.34                  -1.35                -1.28
          N.Valid    777.00                 777.00                 777.00               777.00
        Pct.Valid    100.00                 100.00                 100.00               100.00

Table: Table continues below

 

                    bumbles   clonesize   fruitmass   fruitset   honeybee   MaxOfLowerTRange
----------------- --------- ----------- ----------- ---------- ---------- ------------------
             Mean      0.28       18.77        0.45       0.50       0.42              59.31
          Std.Dev      0.07        7.00        0.04       0.08       0.98               6.65
              Min      0.00       10.00        0.31       0.19       0.00              50.20
               Q1      0.25       12.50        0.42       0.45       0.25              55.80
           Median      0.25       12.50        0.45       0.51       0.25              62.00
               Q3      0.38       25.00        0.48       0.56       0.50              66.00
              Max      0.58       40.00        0.54       0.65      18.43              68.20
              MAD      0.00        0.00        0.04       0.08       0.00               9.19
              IQR      0.13       12.50        0.06       0.11       0.25              10.20
               CV      0.23        0.37        0.09       0.16       2.35               0.11
         Skewness      0.15        0.57       -0.10      -0.52      16.70              -0.02
      SE.Skewness      0.09        0.09        0.09       0.09       0.09               0.09
         Kurtosis      1.68       -0.64       -0.47       0.05     296.31              -1.35
          N.Valid    777.00      777.00      777.00     777.00     777.00             777.00
        Pct.Valid    100.00      100.00      100.00     100.00     100.00             100.00

Table: Table continues below

 

                    MaxOfUpperTRange   MinOfLowerTRange   MinOfUpperTRange    osmia   RainingDays
----------------- ------------------ ------------------ ------------------ -------- -------------
             Mean              82.28              28.69              49.70     0.56         18.31
          Std.Dev               9.19               3.21               5.60     0.17         12.12
              Min              69.70              24.30              39.00     0.00          1.00
               Q1              77.40              27.00              46.80     0.50          3.77
           Median              86.00              30.00              52.00     0.63         16.00
               Q3              89.00              30.00              52.00     0.75         24.00
              Max              94.60              33.00              57.20     0.75         34.00
              MAD              12.75               4.45               7.71     0.18         18.13
              IQR              11.60               3.00               5.20     0.25         20.23
               CV               0.11               0.11               0.11     0.30          0.66
         Skewness              -0.01              -0.01              -0.02    -0.92         -0.21
      SE.Skewness               0.09               0.09               0.09     0.09          0.09
         Kurtosis              -1.35              -1.35              -1.34     0.54         -1.24
          N.Valid             777.00             777.00             777.00   777.00        777.00
        Pct.Valid             100.00             100.00             100.00   100.00        100.00

Table: Table continues below

 

                      Row#    seeds     yield
----------------- -------- -------- ---------
             Mean   388.00    36.12   6012.85
          Std.Dev   224.44     4.38   1356.96
              Min     0.00    22.08   1637.70
               Q1   194.00    33.12   5124.85
           Median   388.00    36.17   6107.38
               Q3   582.00    39.24   7022.19
              Max   776.00    46.59   8969.40
              MAD   287.62     4.56   1417.36
              IQR   388.00     6.12   1897.33
               CV     0.58     0.12      0.23
         Skewness     0.00    -0.06     -0.32
      SE.Skewness     0.09     0.09      0.09
         Kurtosis    -1.20    -0.40     -0.39
          N.Valid   777.00   777.00    777.00
        Pct.Valid   100.00   100.00    100.00

Partición inicial

split <- initial_split(data, prop = 0.70, strata = yield)

Partición inicial

split <- initial_split(data, prop = 0.70, strata = yield)
data_train <- split %>% training()

'data.frame':   541 obs. of  18 variables:
 $ Row#                : int  1 3 4 5 6 12 13 15 16 17 ...
 $ clonesize           : num  37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 ...
 $ honeybee            : num  0.75 0.75 0.75 0.75 0.75 0.25 0.25 0.25 0.25 0.25 ...
 $ bumbles             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ andrena             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ osmia               : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ MaxOfUpperTRange    : num  86 94.6 86 86 94.6 86 94.6 86 86 94.6 ...
 $ MinOfUpperTRange    : num  52 57.2 52 52 57.2 52 57.2 52 52 57.2 ...
 $ AverageOfUpperTRange: num  71.9 79 71.9 71.9 79 71.9 79 71.9 71.9 79 ...
 $ MaxOfLowerTRange    : num  62 68.2 62 62 68.2 62 68.2 62 62 68.2 ...
 $ MinOfLowerTRange    : num  30 33 30 30 33 30 33 30 30 33 ...
 $ AverageOfLowerTRange: num  50.8 55.9 50.8 50.8 55.9 50.8 55.9 50.8 50.8 55.9 ...
 $ RainingDays         : num  1 1 24 34 24 1 16 24 34 24 ...
 $ AverageRainingDays  : num  0.1 0.1 0.39 0.56 0.39 0.1 0.26 0.39 0.56 0.39 ...
 $ fruitset            : num  0.444 0.408 0.354 0.31 0.284 ...
 $ fruitmass           : num  0.425 0.409 0.383 0.366 0.352 ...
 $ seeds               : num  33.4 31.6 28.9 27.3 26.1 ...
 $ yield               : num  4948 4304 3436 2825 2625 ...

Partición inicial

split <- initial_split(data, prop = 0.70, strata = yield)
data_train <- split %>% training()
data_test <- split %>% testing()

'data.frame':   236 obs. of  18 variables:
 $ Row#                : int  0 2 7 8 9 10 11 14 19 27 ...
 $ clonesize           : num  37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 37.5 25 ...
 $ honeybee            : num  0.75 0.75 0.75 0.75 0.75 0.75 0.25 0.25 0.25 0.25 ...
 $ bumbles             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ andrena             : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ osmia               : num  0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 0.25 ...
 $ MaxOfUpperTRange    : num  86 94.6 94.6 77.4 77.4 69.7 86 94.6 77.4 94.6 ...
 $ MinOfUpperTRange    : num  52 57.2 57.2 46.8 46.8 42.1 52 57.2 46.8 57.2 ...
 $ AverageOfUpperTRange: num  71.9 79 79 64.7 64.7 58.2 71.9 79 64.7 79 ...
 $ MaxOfLowerTRange    : num  62 68.2 68.2 55.8 55.8 50.2 62 68.2 55.8 68.2 ...
 $ MinOfLowerTRange    : num  30 33 33 27 27 24.3 30 33 27 33 ...
 $ AverageOfLowerTRange: num  50.8 55.9 55.9 45.8 45.8 41.2 50.8 55.9 45.8 55.9 ...
 $ RainingDays         : num  16 16 34 16 1 16 16 1 16 16 ...
 $ AverageRainingDays  : num  0.26 0.26 0.56 0.26 0.1 0.26 0.26 0.1 0.26 0.26 ...
 $ fruitset            : num  0.411 0.384 0.247 0.428 0.464 ...
 $ fruitmass           : num  0.408 0.399 0.343 0.415 0.436 ...
 $ seeds               : num  31.7 30.5 25 32.3 34.8 ...
 $ yield               : num  3813 3867 2380 4235 5357 ...

Receta de preprocesamiento

data_rec <- recipe(formula = `yield` ~ .,
                   x = data_train) %>%
  step_select(where(base::is.numeric)) %>%
  step_mutate_at(where(base::is.numeric),
                 -contains(c("RainingDays",
                             "AverageRainingDays",
                             "honeybee",
                             "bumbles",
                             "andrena",
                             "osmia",
                             "fruitset",
                             "fruitmass")),
                 fn = ~log(.))
                 
data_prep <- prep(data_rec)
data_train_prep <- juice(data_prep)
data_test_prep  <- bake(data_prep, new_data = data_test)

Definir función para recuperar coeficientes

get_glm_coefs <- function(x) {
  x %>% 
    # get the glm model object
    extract_fit_engine() %>% 
    # transform its format
    broom::tidy()
}

Definir validación cruzada

set.seed(234)
data_folds <- rsample::vfold_cv(data_train,
                                 v = 20)

Regresión múltiple

Definir especificaciones del modelo

# Inicializar un objeto de regresión lineal
linear_model_spec <- linear_reg() %>% 
  # Establecer el modelo de motor
  set_engine('lm') %>% 
  # Establecer el modo de modelo
  set_mode('regression')

Definir fórmula del modelo

formula_lm = yield ~ fruitmass + fruitset + seeds + osmia + andrena + bumbles +
  MinOfUpperTRange + clonesize + AverageRainingDays + honeybee + RainingDays

Definir flujo de trabajo

lm_wf <- workflow() %>%
  add_recipe(data_rec) %>%
  add_model(linear_model_spec,
            formula = formula_lm)

Validación cruzada

tictoc::tic()
validacion_lm <-
  lm_wf %>%
  fit_resamples(
    resamples = data_folds,
    metrics = multi_met,
    control = control_resamples( 
      extract = get_glm_coefs,
      save_pred = TRUE,
      verbose = FALSE
    )
  )
tictoc::toc()

13.86 sec elapsed

Métricas del rendimiento predictivo en fase de validación cruzada

# Métricas promedio de todas las particiones
metrics_lm <- 
  validacion_lm %>%
  unnest(.predictions) %>%
  group_by(id) %>%
  multi_met(truth = `yield`,
            estimate = `.pred`) %>%
  group_by(.metric) %>%
  dplyr::summarise(mean = mean(`.estimate`, na.rm = T),
            std_err = sd(`.estimate`, na.rm = T),
            n = n()) %>%
  ungroup() %>%
  dplyr::mutate(modelo = "lm")
metrics_lm

.metric	mean	std_err	n	modelo
mae	119.4882	15.619637	20	lm
mape	2.1658	0.385811	20	lm
rmse	151.1992	27.588467	20	lm
rsq	0.9881	0.005662	20	lm

Entrenamiento del modelo final

# Métricas promedio de todas las particiones
lm_model <- linear_model_spec %>%
  fit(formula_lm,
      data = data_train_prep)
lm_model

parsnip model object


Call:
stats::lm(formula = yield ~ fruitmass + fruitset + seeds + osmia + 
    andrena + bumbles + MinOfUpperTRange + clonesize + AverageRainingDays + 
    honeybee + RainingDays, data = data)

Coefficients:
       (Intercept)           fruitmass            fruitset               seeds  
          1.984826           -5.944940            2.612865            2.204394  
             osmia             andrena             bumbles    MinOfUpperTRange  
          0.070036            0.063956           -0.058212            0.002545  
         clonesize  AverageRainingDays            honeybee         RainingDays  
          0.031348           -0.035562            0.003223           -0.000328

term	estimate	std.error	statistic	p.value
(Intercept)	1.9848262	0.4838516	4.1021	0.000047388210623994374997604295263187168529839255
fruitmass	-5.9449399	0.3819019	-15.5667	0.000000000000000000000000000000000000000000002951
fruitset	2.6128649	0.1722453	15.1694	0.000000000000000000000000000000000000000000204388
seeds	2.2043935	0.1860721	11.8470	0.000000000000000000000000000069942530762398219600
osmia	0.0700361	0.0100628	6.9599	0.000000000010115665607954346622921262821570564938
andrena	0.0639560	0.0090381	7.0762	0.000000000004724348408114360337636122855187181813
bumbles	-0.0582122	0.0264683	-2.1993	0.028287426461342017602884624238868127577006816864
MinOfUpperTRange	0.0025452	0.0194022	0.1312	0.895680124173958036593035103578586131334304809570
clonesize	0.0313481	0.0065414	4.7923	0.000002145610420601447544935308497926484960771631
AverageRainingDays	-0.0355625	0.0603067	-0.5897	0.555647841162889788435563787061255425214767456055
honeybee	0.0032228	0.0026613	1.2110	0.226440579572288552823167151473171543329954147339
RainingDays	-0.0003279	0.0008122	-0.4037	0.686618002087163170621408880833769217133522033691

Testeo interno

Resumen de la predicción

predicciones <- lm_model %>%
  predict(
    new_data = data_test_prep,
    type = "numeric"
  )
  
predicciones_lm <- predicciones %>% 
  bind_cols(data_test_prep %>% dplyr::select(yield)) %>%
  dplyr::mutate(modelo = "lm",
                yield = exp(yield),
                .pred = exp(.pred))

summary(predicciones_rf)

     .pred          yield         modelo         
 Min.   :2452   Min.   :2380   Length:236        
 1st Qu.:5084   1st Qu.:5118   Class :character  
 Median :6068   Median :6104   Mode  :character  
 Mean   :6061   Mean   :6060                     
 3rd Qu.:7045   3rd Qu.:7035                     
 Max.   :8979   Max.   :8969

Testeo interno

Resumen de las métricas de predicción

error_test_lm  <- multi_met(
  data     = predicciones_lm,
  truth    = yield,
  estimate = .pred,
  na_rm    = TRUE
) %>%
  dplyr::mutate(
    modelo = "lm"
  )
error_test_lm

# A tibble: 4 × 4
  .metric .estimator .estimate modelo
  <chr>   <chr>          <dbl> <chr> 
1 rmse    standard     140.    lm    
2 rsq     standard       0.990 lm    
3 mape    standard       1.89  lm    
4 mae     standard     108.    lm

Testeo interno

Discusión final

Éxito = Conocimiento del Negocio + Conocimiento estadístico + Capacidad de resolver problemas.
Cada modelo necesita de una receta de preprocesamiento adecuada.
Las predicciones deben estar en el mismo dominio que la variable observada.
Muchas variables => maldición de la dimensionalidad => sobreajuste.
No existe una sola forma de resolver un problema.

Aplicaciones de Agro Data Science

Agro Data Science

La agricultura moderna

Agro Data Science

Problemas que resuelven el ADS

R

Tidyverse

Tidymodels

Análisis de mortalidad de insectos

Entendiendo el problema

Partición inicial

Partición inicial

Partición inicial

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Receta de preprocesamiento

Random Forest

Definir especificaciones del modelo

Definir fórmula del modelo

Definir flujo de trabajo

Extraer el conjunto de hiperparámetros del modelo

Definir métricas del performance predictivo

Definir esquema de validación cruzada

Búsqueda del mejor conjunto de Hiperparámtros mediante validación cruzada

Métricas del rendimiento predictivo en fase de validación cruzada

Seleccionar el mejor conjunto de hiperparámetros

Métricas del rendimiento predictivo para el mejor conjunto de hiperparámetros

Entrenamiento del modelo final

Importancia de variables

Testeo interno

Resumen de la predicción

Testeo interno

Resumen de las métricas de predicción

Testeo interno

Predicción de rendimiento

Fórmulas para la proyección del rendimiento

Wild Blueberry Yield Prediction

Partición inicial

Partición inicial

Partición inicial

Receta de preprocesamiento

Definir función para recuperar coeficientes

Definir validación cruzada

Regresión múltiple

Definir especificaciones del modelo

Definir fórmula del modelo

Definir flujo de trabajo

Validación cruzada

Métricas del rendimiento predictivo en fase de validación cruzada

Entrenamiento del modelo final

Testeo interno

Resumen de la predicción

Testeo interno

Resumen de las métricas de predicción

Testeo interno

Discusión final

¿Preguntas?